# This uses the original explicit tau-leaping method by Gillespie -- the
# reason we do not use the Direct Method is not only efficiency, but flexibility.
# Explicit tau-leaping can easily be adjusted to run as the explicit tau-
# leaping method, Euler–Maruyama for real amounts of molecules (CLE), and the
# deterministic Euler method.

model.tau_leaping <- function(a, tau)
	rpois(length(a), a*tau)

model.maruyama <- function(a, dt)
	rnorm(length(a), a*dt, sqrt(a*dt))

model.euler <- function(a, dt)
	a*dt

model.solve <- function(tfinal, X0, V, transitions, params, method, dt)
{
	use_binomial <- identical(method, model.tau_leaping)
	#if(use_binomial) {
	#	binomial <- apply(V<0, 1, any)
	#	Vleft <- V<0
	#}
	
	Xt <- matrix(c(time=0,X0), 1,length(X0)+1)
	t <- 0
	X <- X0
	while(t < tfinal) {
 		a <- transitions(X)
#  		for(r in 1:nrow(V)) {
#  			v <- V[r,]
#  			if(use_binomial & any(v < 0)) {
#  				n <- sum((v<0)*X)
#  				lambda <- a[r]
#  				if(n > 0)
#  					k <- rbinom(1, n, 1-exp(-lambda*dt/n))
#  				else
#  					k <- 0
#  				#k <- rbinom(1, n, 0.5)
#  			}
#  			else
#  				k <- method(a[r], dt)
#  			X <- X + as.vector(k*v)
#  		}
		
		k <- method(transitions(X), dt)
		X <- X + as.vector(k %*% V)
		X[X<0] <- 0  # dont allow negative

		t <- t + dt  # dt is aka tau
		Xt <- rbind(Xt, c(time=t, X))
	}
	Xt
}
